#ifndef _GRID_H_
#define	_GRID_H_

#include <vector>
#include <algorithm>
#include <functional>
#include <cmath>

/*
** Low-level rectangular grid to store value at grid vertices
** ValueType: the data type the vertices store
*/
namespace Packing_common
{

	struct Local_frame_2
	{
		typedef CGAL::Vector_2<CGAL::Exact_predicates_inexact_constructions_kernel>	Vector_2;
		typedef CGAL::Point_2<CGAL::Exact_predicates_inexact_constructions_kernel>	Point_2;
		typedef CGAL::Aff_transformation_2<CGAL::Exact_predicates_inexact_constructions_kernel>	Transformation_2;

		// x_base, y_base, origin are base vectors and origin of this frame in **world coordinate**
		Local_frame_2(const Vector_2& x_base, const Vector_2& y_base, const Point_2& origin) : u(x_base), v(y_base), o(origin) 
		{
			u = u / CGAL::sqrt(u.squared_length());
			v = v / CGAL::sqrt(v.squared_length());
		}
		//Local_frame_2()	{ }
		
		inline void transform(const Transformation_2& t)
		{
			u = t(u);
			v = t(v);
			o = t(o);
		}

		inline Point_2 local_coord(const Point_2& world_coord)	const
		{
			Vector_2 op(o, world_coord);
			return Point_2(op*u, op*v);
		}

		inline Point_2 world_coord(const Point_2& local_coord) const
		{
			return o + local_coord.x()*u + local_coord.y()*v;
		}
		inline Vector_2 world_coord(const Vector_2& local_coord) const
		{
			return local_coord.x()*u + local_coord.y()*v;
		}

		Vector_2 u;
		Vector_2 v;
		Point_2 o;
	};

	template <class ValueType>
	class Grid_2
	{
	public:
		Grid_2(double x_min, double x_max, double y_min, double y_max, size_t nrows, size_t ncols) ;
		~Grid_2();

		//access
		ValueType& at(size_t row_id, size_t col_id) ; // exact query to vertex value by grid index

		//template <class InterpolationPolicy>
		//ValueType at(double x_coord, double y_coord, InterpolationPolicy interp_policy) ; // approximate query to any point inside the grid by interpolation

		const ValueType& at(size_t row_id, size_t col_id) const;

		template <class Operator>
		void for_each_node(Operator op);

		// query
		inline size_t nb_rows() const { return nrows; }
		inline size_t nb_columns() const { return ncols; }
		inline double xmin() const { return min_x; }
		inline double xmax() const { return max_x; }
		inline double ymin() const { return min_y; }
		inline double ymax() const { return max_y; }
		inline double xstep() const { return x_step; }
		inline double ystep() const { return y_step; }

	private:
		std::vector< std::vector<ValueType> > values; // the values stored at every vertex
		size_t nrows;
		size_t ncols;
		double min_x;
		double max_x;
		double min_y;
		double max_y;
		double x_step;
		double y_step;
	};

	template <class ValueType>
	Grid_2<ValueType>::Grid_2(double x_min, double x_max, double y_min, double y_max, size_t row_cnt, size_t col_cnt)
	{
		assert(x_min < x_max && y_min < y_max);
		min_x = x_min;
		max_x = x_max;
		min_y = y_min;
		max_y = y_max;
		nrows = row_cnt;
		ncols = col_cnt;
		x_step = (max_x - min_x)/ncols;
		y_step = (max_y - min_y)/nrows;
		values.resize(nrows+1);
		std::for_each(values.begin(), values.end(), std::bind2nd(std::mem_fun_ref(&std::vector<ValueType>::resize), ncols+1));
	}
	template <class ValueType>
	Grid_2<ValueType>::~Grid_2()
	{
		std::for_each(values.begin(), values.end(), std::mem_fun_ref(&std::vector<ValueType>::clear));
		values.clear();
	}

	template <class ValueType>
	ValueType& Grid_2<ValueType>::at(size_t row_id, size_t col_id)
	{
		return values.at(row_id).at(col_id);
	}

	template <class ValueType>
	const ValueType& Grid_2<ValueType>::at(size_t row_id, size_t col_id) const
	{
		return values.at(row_id).at(col_id);
	}

	template <class ValueType> template <class Operator>
	void Grid_2<ValueType>::for_each_node(Operator op)
	{
		for (size_t i = 0; i < values.size(); i++)
			std::for_each(values[i].begin(), values[i].end(), op);
	}
	
}
#endif